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ABSTRACT 

The Kepler mission has discovered about a dozen circumbinary planetary systems, all containing 
planets on ^ 1 AU orbits. We place bounds on the locations in the circumbinary protoplanet ary disk, 
where these planets could have formed through collisional agglomeration starting from small (km-sized 
or less) planetesimals. We first present a model of secular planetesimal dynamics that accounts for 
the (1) perturbation due to the eccentric precessing binary, as well as the (2) gravity and (3) gas drag 
from a precessing eccentric disk. Their simultaneous action leads to rich dynamics, with (multiple) 
secular resonances emerging in the disk. We derive analytic results for size-dependent planetesimal 
eccentricity, and demonstrate the key role of the disk gravity for circumbinary dynamics. We then 
combine these results with a simple model for collisional outcomes and find that in systems like Kepler 
16, planetesimal growth starting with 10-100 m planetesimals is possible outside a few AU. The exact 
location exterior to which this happens is sensitive to disk eccentricity, density and precession rate, as 
well as to the size of the first generation of planetesimals. Strong perturbations from the binary in the 
inner part of the disk, combined with a secular resonance at a few AU inhibit the growth of km-sized 
planetesimals within 2 — 4 AU of the binary. In situ planetesimal growth in the Kepler circumbinary 
systems is possible only starting from large (few-km-sized) bodies in a low-mass disk with surface 
density < 500 g cm“^ at 1 AU. 

Subject headings: planets and satellites: formation — protoplanetary disks — planetary systems — 
binaries: close — Kepler 16 — accretion disks 


1. INTRODUCTION 

Recent discoveries of exoplanets in stellar binaries by 
radial velocity surveys and the Kepler mission stimulated 
significant interest in understanding the origin of such 
planetary systems. Planet-hosting stellar binaries come 
in two flavors , commonly denoted as S-type or P-type 
(Dvorak 1982|). They correspond to systems where the 
planet orbits one star with the other as an external per- 
turber (S-type), or where the planet is in orbit around 
both stars (P-type). 

In this paper, we focus on P-type or circumbinary sys¬ 
tems. Kepler has so far revealed to us eight such systems, 
all containing sub-Jovian planets in orbits around main- 
sequence binaries. These systems have a range of param¬ 
eters, but generally the stars have masses ^ Mq, and 
are on moderately eccentric orbits with semi-major axes 
ab ~ 0.1 — 0.2 AU. The planetary orbits have semi-major 
axes smaller than 1.1 AU. The parameters of known Ke¬ 
pler circumbinary planets are summarized in Table 

There may also exist another population of circumbi¬ 
nary planets hinted at by transit timing of post-common 
envelope binaries. The most plausible of these are two 
planets arou nd the NN Serpentis binar y system (Marsh 


et al . 20141. It has been suggested (Volschow et al. 


2014p that such planets are not primordial and formed 


from matter ejected during common envelope evolution. 
There are also a few directly imaged long period plan¬ 
etary mass circumbinary companions at projected sepa- 
rations > 100 AU from the host star, e.g. ROXs 42Bb 
(Currie et al. 2014). We will not address the origin of 
such systems in this work. 
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Planet formation in a circumbinary disk faces chal¬ 
lenges due to vigorous planetesimal excitation driven by 
the non-Newtonian potential of the binary. Large rela¬ 
tive velocities between planetesimals are harmful to co¬ 
agulation because collisions lead to planetesimal destruc¬ 
tion rather than merging. It is generally believed that the 
circumbinary planets cou ld not have formed in situ via 


et al. 

2012 

speed 

S of K 


(2014) proposes a model for in situ formation, in which 


large planetesimals form at the pressure maximum near 
the inner edge of the disk where t he surface density of 
solids is expected to be increased. Bromley fc Kenyon] 
(2015) suggest that circumbinary planetesimals may set¬ 
tle onto special class of orbits where their growth would 
be promoted, and we will comment on this work further 
(^). 


Recently Rafikov (2013) proposed that the gravita¬ 
tional effect of a rnassive axisymmetric protoplanetary 
disk may strongly suppress eccentricities of circumbinary 
planetesimals. That happens because disk gravity drives 
rapid relative precession of planetesimal and binary or¬ 
bits, reducing the time-averaged non-axisymmetric com¬ 
ponent of the binary potential, which drives planetes¬ 
imal eccentricity. However, simulations show that cir- 


develo 

p eccentricity themselves ( 

Pelupessy & Portegies 

Zwart 

2013 

Meschiari 

2014 

). Silsbee & Rafikov (2015; 


gravitational effect of an eccentric disk, and applied those 
results to planet formation in S-type systems. Rafikov & 
Silsbee (2015a; hereafter RS15a), additionally included 
the effects of gas drag, which damps free eccentricity 
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TABLE 1 

Binary system parameters 


devoid of gas (Artymowicz fc Lubow||1994 Pelupessy & 


Systerr 


M, 


Rs 


db 




K34 

1.05 

1.16 

1.02 

1.09 

0.23 

0.52 

1.09 

K16 

0.69 

0.65 

0.20 

0.23 

0.22 

0.16 

0.70 

K47 

1.04 

0.96 

0.36 

0.35 

0.084 

0.024 

0.30, 0.9 

K38 

0.95 

1.76 

0.25 

0.23 

0.15 

0.10 

0.48 

KIC 4862625 

1.47 

1.7 

0.37 

0.34 

0.18 

0.20 

0.64 

K 413 

0.82 

0.78 

0.54 

0.48 

0.10 

0.037 

0.36 

K35 

0.89 

1.03 

0.81 

0.79 

0.18 

0.14 

0.60 

KIC 9632895 

0.93 

0.83 

0.19 

0.2 

0.18 

0.05 

0.79 


Portegies Zwart 2U13). According to these calculations 
±he inner edge oi the disk is truncated at the inner radius 
Oin = (1-7 — 3.3)a{,, increasing as Cb increases from 0 to 
0.7. Most of the observed stellar orbits in circumbinary 
systems are on the low end of that eccentricity range 
(see Table [^, so we assume that the disk is truncated 
on the inside at Oin = 2ah. This is somewhat lower than 
Oin « 306 favored by jPelupessy fc Portegies Zwart ( 2013 ). 

Our disk model is analogous to that adopted in |Silsbee| 
& Rafikov (2015). We assume that circumbinary disk 


‘‘References: lOrosz et al.lll2012all:IWelsh et al. 
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uIDovle et al. 

f( 

12011 
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^Mp (Ms) is the mass of primary (secondary) star, Rp (Rs) is 
the radius of primary (secondary) star, a;, is the binary semi-major 
axis, Cb is the binary eccentricity, ap is the planetary semi-major 
axis. 

‘^This system contains at least two planets. 


treamlmes are ellipses with a common apsidal line and 


which may be a function of time (see Section 5.1) if the 
disk precesses as a solid body (see a discussion of this 
assumption in l|^. We adopt power law scalings for the 
surface density at periastron and eccentricity of gas 
streamlines given by 


and produces apsidal alignment of equal-size planetes- 
imals. The dynamical results of RS15a were applied in 
Rafikov and Silsbee (2015b; hereafter RS15b) to study 
collision outcomes in S-type planetary systems. 

In this paper we extend these calculations to P-type 
systems. In the following discussion of planetesimal dy¬ 
namics we account for the combined effects of (1) gas 
drag, (2) gravitational perturbations from the central bi¬ 
nary and (3) gravitational perturbations from the eccen¬ 
tric protoplanetary disk, and derive expressions for the 
encounter velocities of circumbinary planetesimals. We 
then use these results to determine which sizes of plan¬ 
etesimals are able to grow at different locations in the 
disk through mutual collisions. 

In Section!^ we describe our model system. In Section 

we write down the disturbing function characterizing 
the gravitational perturbations from the binary and the 
disk, and present the equations governing planetesimal 
dynamics. In Sections i-[zl we describe planetesimal 
dynamics in several different regimes. Section [^discusses 
collision outcomes calculated using the prescription de¬ 
scribed in Appendix Section P critically assesses our 
underlying assumptions, and explores the outcome of re¬ 
laxing some of them. Section compares our results 
with those in the literature. Finally, we summarize our 
main conclusions about circumbinary planet formation 
in Section El 


2. GENERAL SETUP 

Our model system is a close stellar binary consisting 
of a primary with mass Mp and a secondary with mass 
Ms < Mp in orbit with semi-major axis at and eccentric¬ 
ity Cfc. Orientation of the binary is given by the apsidal 
angle Wb with respect to a fixed reference direction. Bi¬ 
nary orientation (and Wb) may in general vary in time as 
a result of binary precession; we consider this possibility 
in Section!^ For convenience, we define Mb = Mg -b Mp, 
and /r = Mg/Mb- Throughout this paper, unless other¬ 
wise noted, we provide numerical estimates for a fiducial 
system which has the binary parameters of Kepler 16: 
Mp = 0.69 Mq, Mg = 0.2 Mq, at = 0.22 AU, = 0.16. 

Orbiting the barycenter of the binary is a gaseous pro- 
toplanetary disk of mass Md- Simulations generally find 
disks around binaries with fj, ^ 1/2 to be tidally trun¬ 
cated on the inside, resulting in the inner cavity relatively 


Sd(a) = So , ed(a) = eg 


_ /ooy 


( 1 ) 


where a is the semi-major axis of the elliptical fluid tra¬ 
jectory, and Oq is a reference distance, which we take to 
be 1 AU. These dependencies lead to non-tr ivial surface 
density behavior described in Statler (2001); SR15. 


TABLE 2 

Fiducial system parameters 


Parameter 

Value 

Mp 

O.69M0 

Ms 

O.2OM0 

P 

0.22 

db 

0.22 AU 

ef, 

0.16 

So 

3,000 g cm“2 

eo 

0.024 

p 

1.5 

9 

1 


For numerical estimates in this work we adopt the 
following set of disk parameters: Sq = 3,000 g cm“^, 
eo = 0.024, and and Cd slopes p = 1.5 and <7=1. 
These and Kepler-16 binary parameters are listed in Ta¬ 
ble [H for convenience. 

Our choice o f p corresponds to the MMSN model of 
Hayashi (1981). This value of p also follows from the de- 
cretion disk iriodel of Pringle (1991), assuming no mass 
to pa ss through the inner bouiidary of the disk (Rafikov 
2013). Our fiducial values of q and eg are chosen to 


correspond to the behavior of the forced eccentricity of 
free particles or biting in the potential of the binary with 
eccentricity (Moriwaki & Nakagawa 2004), see Equa¬ 
tion (10). This estimate is just a reasonable zeroth order 
guess for ed{a), as eccentricity of the fluid disk is also 
affected by pressure, viscous forces and disk gravity. 

Assuming the mass to be concentrated in the outer 
part of the disk (i.e. p < 2), we can relate Md to the 
outer radius of the disk Oout and Eg: 


27r 

Md=- -Egoga; 


2-p 

i 0.037 Mq 


out 


-'0 / ^out \ 


0.5 


3,000 gcm-2 V75AU/ 


( 2 ) 

( 3 ) 
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where the numerical estimate has been performed for 
p = 1.5. The outer truncation radius of 75 AU is char- 
acteristic of the more massive disks aro und single stars 
( Andrews et al.|[2009 Harris et al.||2012 1. 

3. EQUATIONS OF PLANETESIMAL DYNAMICS 

We now outline the mathematical framework that we 
use to describe planetesimal dynamics. A planetesimal 
orbit is characterized by the semi-major axis Cp with re¬ 
spect to the binary barycenter, eccentricity Cp, and apsi- 

dal angle vjp; its mean motion is rip = y^GM^/o^. 

Planetesimal motion is affected by both conservative 
forces — gravity of the binary and the disk — a n d non - 
conservative gas drag. Analogous to Rafikov (2013), 
SR15, RSlSa, we employ secular perturbation theory 
(Murray & Dermott||1999 Moriwaki & Nakagawa||2004 1 
to determine the effect ot the former on planetesiinal ec- 
centricity behavior. This approach uses the disturbing 
function, described next in §3.1[ which accounts for the 
perturbations to planetesimal inotion produced by the 
gravity of a massive protoplanetary disk and the non- 
Newtonian gravity of the binary. Our treatment of the 
effe cts o f gas drag on planetesimal dynamics is d escri bed 
in §3.2[ All these components are combined in §3.3[ re¬ 
sulting in a set of general equations (18l-(19) describing 


secular planetesimal dynamics in circumbinary planetes¬ 
imal disks. 

3.1. Disturbing Function 

Planetesimal disturbing function R — Rd + Rb consists 
of contributions due to an eccentric disk Rd and due to 
the non-Newtonian gravity of the binary Rb- 

SR15 derived the disturbing function due to a disk with 
eccentricity and surface density given by Equation ([^, 
and distance-independent apsidal angle Wd as 


Rd — npQj 


.pu,p 


Ad = 2'K'ipi 


^Cp -I- BdCp cos (tUp - Wd) 

GY.d{ap) 

QjpTip 

■yr 

So V'l 


(4) 

(5) 


!-1.6 X 


Bd = T^i’2 


3,000 g cm-2 (-0.55) 

G'Fid{Op')Cd{p,p') 


0.89Mo 

Mb 


0.5 


P.5’ 


( 6 ) 


LlpILp 

1.3 X 10"®yr"^ 


'ip2 ( O.S9 Mq 


3,000 g cm-2 1 85 


Mb 


0.5 


eo -1 


0.024 


and Op^s = ap/(5AU). Here ipi and '02 are coefficients of 
order unity (SR15), which are effectively constants far 
from the disk edges, for Oin ^ cip < Cout- For p = 1.5 
and q = 1, SR15 show that '0i « —0.55 and tp 2 ~ 1-85, 
except near the edges of the disk. As shown in SR15, Ad 
and Bd are dominated by local disk properties for our 
choice of p and q, so non-power-law behavior of surface 
density or eccentricity near the edges of the disk will not 
greatly affect th e values of •0i and 02- 


According to Moriwaki & Nakagawa (2004) the dis¬ 
turbing function due to the binary has the form similar 


to Equation (|^: 


Rb = Una. 


p“p 


A 


5 „2 


2 + BbGp cos {Wp - Wb) 


4 Up \ Op 


(7) 

( 8 ) 


Mb 


0.89Mq 


0.5 


a, 


0.17 V0.22AU7 

,-3.5 
'P.5 > 




i6.5 X 10"V"^ 


Ob 


(9) 


(a 


.V 


Mb 


O.89M0 


0.5 


0.096 V0.22AU7 

-6^ -4.5 
0.16 ■ 


Here Ub = \jGMb/a^ is the mean motion of the binary, 
and /(^) = pl{ 1 — p){l — 2/rb We use these disturb¬ 
ing functions in Sections |^ - |^ to calculate the orbital 
dynamics of planetesimals in several different dynamical 
regimes. 

In the absence of gas drag and disk gravity, free par- 
ticles in the binary potential attain forced eccentricity 
( Moriwaki fc Nakagawa|[2004 | 


__ 5 , .Ob 

^forced — ^ 

4 ap 


( 10 ) 


! 0.024 


ab 


lAU 


(1 - 2/4) Cb 
0.56 0.16 0.22AU 


which provides a useful reference value, e.g. for our 
choice of eg. 

Introducing planetesimal eccentricity vector Cp = 
(fcp, hp) = ep{cos Wp, sin Wp) we can then write down the 
full disturbing function as 


R = Upttp 


^ (/ip + fcp) + Bdkp cos xud 


2 v^p ' ""pJ ' -^a-'^p 

-I- Bdhp sin Wd + Bbkp cos Wb + Bbhp sin Wb] 


( 11 ) 


where Wd = w^dit), Wb = Wb(t). 

Here A = Ad + Ab is a. precession rate of planetesimal 
free eccentricity due to the axisymmetric components of 
both the disk and binary gravity. Its behavior as a func¬ 
tion of distance from the binary is shown in Figure to¬ 
gether with separate curves for Ad{ap) and Ab(ap). The 
fact that A goes through zero at some semi-major axis 
a A has very impor tant implications for planetesimal dy¬ 
namics, see §3.1.1| 


The other terms in equation (11) that depend on the 


planetesimal orientation (i.e. Wp) describe excitation of 
planetesimal eccentricity by the torques produced by the 
non-axisymmetric components of the disk and binary po¬ 
tentials. We discuss the relative role of different contri¬ 
butions to R next. 

3.1.1. Transition Between Binary-Dominated and 
Disk-Dominated Regimes 
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Fig. 1. — Planetesimal precession rate A = Afj A^ due to both 
binary and disk gravity (black curve) as a function of calcu¬ 
lated assuming the fiducial system parameters in Table Dot¬ 
ted and dashed curves represent A(i{ap) and ^^(ap). Several ref¬ 
erence values corre spon ding to disk and binary precession rates 
considered in ^5.1| and ^6.1| are also shown (expressed in units of 
^leAU < 0 — planetesimal precession rate at 6 AU). The vertical 
dashed line marks ap = see equation l |12[ |. For our fiducial 
binary -h disk model (see Table A\qau = —1-2 x 10“^yr“^, 
= 2.6 X 10“^ yr~^- 


Because and Bt fall off more rapidly with semi¬ 
major axis than their disk-related counterparts and 
Bf,, see Equations (|^-(|^ and we find that out¬ 

side a certain radius the disturbing function R should be 
dominated by the disk gravity. This situation is analo¬ 
gous to the so-called DD regime discussed in SR15, in 
which both the axisymmetric and the non-axisymmetric 
components of the disturbing function are dominated by 
the disk. 

In the opposite limit, very close to the binary, the 
gravitational perturbations are dominated by the binary 
system. This is analogous to the so-called BB regime 
of SR15 in which both the axisymmetric and the non- 
axisymmetric components of the disturbing function are 
dominated by the binary. At the same time, the disk still 
cannot be ignored because of the gas drag. 

Equations ® and ® predict that planetesimal pre¬ 
cession (i.e. switaies from being binary- to disk- 
dominated at a characteristic semi-major axis a a, where 
|A,| = |A,|: 


a A = 


3/r(l - fi) M^al 
87r|'0i| 

!l.9AU ' 


l/(4-p) 


( 12 ) 


(a 


V 


Mb 


22AU/ 0.89Mo 

3,000gcm“^ / /x(l — 
I 0.17 


2/5 


where the numerical estimates have been performed for 
our fiducial system. This transition occurs via a secular 
resonance where A = 0 , emerging because disk and bi¬ 
nary drive planetesimalprecession in opposite directions, 
see equations ([^ and 



Fig. 2.— Three different dynamical regimes in the space of disk 
density and semi-major axis. Calculation is done for our fiducial 
system parameters {Mp = 0.69Mq, Ms = 0.2Mq, aj, = 0.22 AU, 
= 0.16, p = 1.5, and q = 1), but with lowered disk eccentcity 
at 1 AU eo = 2.4 x 10“^, which is 10% of the forced eccentricity 
to broaden the DB regime. The vertical red line shows the 
location at which binary precession rate equals to the disk-driven 
plan etesimal precession, the significance of which is discussed in 

3631 


On the other hand, planetesimal eccentricity excitation 
switches from being binary- to disk-dominated at the dif¬ 
ferent characteristic distance as, where \Bd\ = \Bb\'. 


as = 


15/(/t) Mbal Cb 
16TT\'ip2 \ eo 

1.5AU 


1 l/(5-p-ij) 


(13) 


eb 

0.024 

Mb 


0.16 

eo 

0.89Mq 

V 0 . 22 AUy 


3000 g cm ^ /(/i) 
Eq 0.096 


2/5 


For a single component of the disturbing function (ax¬ 
isymmetric or non-axisymmetric), the region of transi¬ 
tion between binary-domination and disk-domination is 
quite narrow, since both A^/Ab and Bd/Bb are rising fast 

with Op (as Op^^). 

If a A is widely separated from ub, there will also be a 
region of the disk where excitation is dominated by the 
disk and precessio n by the binary, or vice-versa. In par¬ 
ticular, in Section [ 6 .2[ we consider the limit of a low ec¬ 
centricity disk (eo —> 0 ) in w hich qa c an b e substantially 
less than os, see Equations (12), and (13). Planetesimal 
dynamics in the intermediate region qa cip as are 
then analogous to the DB regime of SRI 5. In the oppo¬ 
site limit of a high-eccentricity disk, eo = ed(l AU) >0.1 
(which is probably not very realistic), there would exist a 
region as Op < 0 ^ 4 , in which the eccentricity excitation 
is dominated by the disk, while precession is controlled 
mainly by the binary — analogous to the BD regime of 
SR15. 


The locations of these regimes in the space of disk den¬ 
sity and eccentricity are illustrated in Figure The 
width of the DB regime depends on the degree to which 
the eccentricity of the disk falls below the free particle 
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eccentricity in the binary potential. Figure also in¬ 
cludes the line where \wb\ = |^d|- This gives roughly the 
radius outside of which excitation due to the bina^ is 
substantially reduced by its precession (see Section]^. 

3.2. Gas drag 

Following RS15a, we consider gas drag to damp plan- 
etesimal eccentricity as 


dep 

dt 


rd 


(14) 


ate (Weidenschilling 
coeflicient, pd is the 


where = (fed, hd) = Cdicos Wd, sin Wd) is the eccentric¬ 
ity vector of the local disk fluid element, and is the 
eccentricity damping time. This implies that gas drag 
drives planetesimal orbits towards full alignment with 
the gas trajectories on a characteristic timescale Td- 
This timescale is not independent of ep in general. 
For planetesimals with radii dp on the order of km, a 
quadratic drag law F = — ((7£)/2)7r(ipPd'Cr-Vr is appropri- 
1977| ), where Cd ~ 0 .5 is the drag 
ocal gas density and is the rela¬ 
tive planetesimal-gas speed. In this case Td oc and 
RSI5a had Td oc where e^. = ep — is the relative 
eccentricity between the object and the gas. 

In this work we have chosen to also account for the fact 
that gas orbits the binary more slowly than a planetesi¬ 
mal because of the radial pressure support in a gaseous 
disk. This gives rise to additional (predominantly az¬ 
imuthal) irreducible relative velocity between the gas 
and planetesimal Av^ = tjvk, where vk is the Keple- 
rian speed and the explicit expression for 77 ^ 1 is given 
by eq. (13) of RS15b. This velocity differential does not 
vanish even when fluid and planetesimal orbits coincide 
and Bp = ed- To describe this effect, we introduce fidu¬ 
cial eccentricity 


— 


Avd 


2E{y/3/2) VK 
,0.0038^ay^ 

Mb 


(15) 


The numerical estimate uses the prescription for the scale 
height h given in RS15b (their eq. (14)), adapted for a 
disk temperatur^ of 200 K at 1 AU: 


= 0.028 


-^0 / Qp \ 

Mb VAU/ 


(16) 


We generalize the expression for the characteristic 
damping time Td from RS15a as follows: 


Td = 


25/2^3/2 

3C'BE(y3/2)'“^ 
1.6 X 10^ yr 


—1 Pp^p ^ 


v- 1/2 


Yjd ap 

Pp dp 


(17) 


3 g cm 3 km 
0.89Mo 3,000 gcm-2 0.01 


Mb 


13/4 
>,5 ■ 


^ Stars in P-type binaries typically have lower masses and lower 
luminosities that the members of S-type binaries studied in SR15b. 


Here, E(-\/3/2) « 1.21 is a complete elliptic integral, pp 
is the planetesimal bulk density, and we have assumed 
p = 1.5 in the estimate. 

Note that the numerical coefficient in equation 0 is 
different from the analogous expression in RS15a because 
they adopted a rough estimate pd = Ed/h, whereas here 
we use Pg = Tid/'/^h, appropriate for an isothermal 
disk. We have picked the coefficient in the definition of 
ed, to give the corre ct damping time (in agreement with 
Adachi et al.| (|1976|)) in the limit that ^ e^. 

Fquation (17| is only approximate for ~ e,.. How¬ 
ever, it does capture the reduction of Td for particles with 
low values of due to the irreducible velocity differential 
Au 0 caused by pressure support. It is thus an improve¬ 
ment over the approximation used in SR15a. Low values 
of Or ^ e^j) are more typical at several AU around P-type 
binaries than in S-type systems. 


3.3. 


General evolution equations 

and |3.2| 


We now combine the results of 13.1.1 


We 


use the standard Lagran ge equations ([Murray fc Der- 


mott 

1999 

Rafikov||2013 

to relate dkp/dt and dhp/dt to 

the d 
due t 
est 0 

isturbmg function 

0 gas drag given by 
rder in eccentricity 

11 

ec 

tt 

). Addi 
nation ( 
le follow 

14 

an, 

the contributions 
), we find, to low- 
1 set of evolution 


equation for planetesimal eccentricity e. 


p- 


dk 

=-Ahp - BdSinwd{t) - Bb sinwbit) 

kp kd 


Td 

= Akp + Bd cos zudit) + Bb cos Wb{t) 
hp hd 


Td 


(18) 


(19) 


As before, vjd and Wd are in general functions of time. 

These master equations provide a basis for subsequent 
analysis of planetesimal dynamics in Sections i-m 

3 . 4 . General Remarks on Planetesimal Dynamics 

Before embarking on a detailed discussion of planetes¬ 
imal dynamics in the following sections, we outline some 
general features of planetesimal eccentricity evolution de¬ 
scribed by equations (18)-(19|. 

First, in the case ^quadratic drag these equations 
do not in general admit an analytical solution for arbi¬ 
trary time dependence of Wd{t) and vubit). At the same 
time, there are several important limits where analytical 
treat ment is possible, and these situations are covered in 
0 |6.2| These solutions allow us to gain important 

insists into how the binary or disk precession may af¬ 
fect planetesimal dynamics, which remain valid in more 
complicated setups (Q. In addition, some features of 
the general planetesirnal dynamics with both Wd and Wb 
varying in time can be gleaned by consideringa simpler 
case of linear gas drag, covered in Appendix and dis¬ 
cussed in 0 

Second, we find quite generally that any free eccen¬ 
tricity describing the initial conditions for planetesimal 
evolution damps away on a characteristic time ~ r^. Asa 
result, planetesimal eccentricity inevitably converges 
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to its forced value, which is determined by many factors, 
see 0. This fe ature of the evolution has been previously 
pointed out in Beauge et al. (2010) and RSlSa, and im¬ 
plies that after the initial transient lasting for ^ plan- 
etesimals lose memory of their initial conditions. 

In particular, collisions between planetesimals, which 
perturb them away from the equilibrium eccentricities, 
may be considered as a minor effect for the dynamics as 
long as they are infrequent enough, i.e. the mean time 
between them is longer than We p rovide a discussion 

Convergence of 


ifies our analysis as 


of this approximation in Section 9.3 
to a certain fixed state greatly simp 
we see in the following sections. 

Circumbinary systems exhibit a wide range of plan- 
etesimal dynamical behavior throughout the disk. We 
now describe them under different assumptions about 
the disk and binary precession. We start by consider¬ 
ing a simple case where we ignore the effects of binary 
precession. Although the latter is likely very important 
in reality, ignoring it at first allows us to better illustrate 
certain aspcts of planetesimal dynamics. Thus, we con¬ 
sider the case of non-precessing disk and binary in ^ 
We then include the possibility of the disk precession in 
^ Finally, in and 0 we consider a possibility of the 
binary precession. 


4. 


PLANETESIMAL DYNAMICS IN THE ABSENCE OF 
DISK AND BINARY PRECESSION 


If we ignore precession of both the disk and the binary, 
i.e. take zud(t) and Wh(t ) to be constant in time, we 
can easily solve equations (|18|)-( 19). This approximation 


should be valid for small planetesimals, which have stop¬ 
ping times shorter than the two precession timescales, 
Td '^di'^b- In this case the eccentricity vector 
should rapidly adjust to the fixed values corresponding 
by the instantaneous values of vud and Wb- 

Solutions for planetesimal dynamics in this limit ac¬ 
counting for both the disk and the binary gravity have 
been previously derived in RS15a, but for S-type bina¬ 
ries. They remain fully valid in the circumbinary case 
as well, as long as At and Bb are understood as being 
represented by our Equations ([^ and (1^. The forced 
eccentricity to which ep converges after the initial period 
of damping the free eccentricity, is given by a sum of 
forced terms due to the binary and to the disk ef^d'- 


= Bf^b + e/,d, 





e/,d = 


cos {rub + 4>b) 
sin {wb + (t>b) 

cos {Wd + (t>d) 
sin(rod -I- (j)d) 


_1 -I- {ATd)^_ 
where (j)d and (pb are phase angles given by 
ed - ABdTj 


cos pd = 


and 


, -ATd 

COS 0b = 


( 20 ) 

( 21 ) 

( 22 ) 

(23) 

(24) 


y/l+ {ATdY 

One can see that in general Bp is misaligned with both 


the disk and the binary as it is parallel to neither nor 
Gfe = eb{cosvub,sinzub)- It is also clear that Bp does not 
vary in time when Wd and vjb are fixed. 


4.1. Relative Planetesimal-Gas Eccentricity 

Solution (31) implies that the relative planetesimal-gas 
eccentricity = |ep — e^il is given by (RS15a) 


Avd 


Ct> — 6c 


(25) 


IS 


0 1 + i^Tdf 

where the characteristic eccentricity e 
Sc = 1^1 ^ [(4led -l- Bd)^ B^ 

-1-2 cos {wd - Wb)Bb{Aed + Bd)f^^ ■ (26) 
We find it convenient to define a characteristic size d, 
for which Ard = 1 when -I- is replaced with Cc ii 
Equation i.e. eccentricity damping time is of order 


the orbit precession time scale. Equation ( |17[ ) implies 
that 


dr = 


3C£iE(-\/3/2) T,d Up Op 


= 1.3m 


25/27r3/2 

Mb 


0.89Mq 


Pp \A\ h 

1.5 


3 g cm i3/4 


Pp 




(27) 


Because the numerical estimate is made at 5 AU, far 
outside of a a and ob for typical disk parameters, we have 
assumed for simplicity that the contribution of the binary 
to the disturbing function is negligible, i.e. A « A^, 
Bb « 0. In this regime (far from the star) both Cc and 
dc are independent of disk mass Md (or Eq). This can 
be seen from Equations ( [26| and ( [27| ), and the fact that 
Ad, Bd and E^ are all proportional to Md- 
Using our newly defined dc, we may rewrite equation 

(17) as Ard = {dp/dg) {ed -I- . Combining this 

result with Equation (|32| we can solve for Ard and e^: 


AdTd = 2 


1-K + L 


1 1/2 


i^dg/dp 


K 


1-K + L 


1 1/2 


1 K -\- 2(dgldp)^ L 


(28) 


(29) 


where 


K = {e^/egf{dg/dpf, 

L=[{K +lf +A{dgldpf]^'^ . 

For most of the paper, we will be considering situations 
in which Cg e^. In this limit, K ^ 0, L ^ 2dc/dp, 
and planetesimal dynamics are determined purely by the 
values of dg and eg, so that our results reduce to those of 
RSI5a, with Td longer by a factor of -x/Stt (to compensate 
for the different definition of pd)- 


4.2. Solution far from the Binary 

A useful limit of planetesimal dynamics without disk 
or binary precession is obtained when we are justified in 
neglecting the binary perturbation. This regime is nat¬ 
urally realized far from the binary. This is equivalent to 
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ap, AU 


dp, AU 


Fig. 3. — Radial dependence of (see Equation | |25[ |) in a sys¬ 
tem with no binary or disk precession, but including gravity of 
both the binary and the disk, for 4 different disk models. Calcu¬ 
lations assume the binary parameters of Kepler 16 with fiducial 
disk parameters unless otherwise labelled; z u^ = = zui, is assumed. 
Characteristic ec is determined by Equation (|26| >. Different colors 
correspond to different planetesimal sizes, as indicated in panel D. 

the DD regime of SRI 5 in which disk gravity dominates 
dynamics, so \Bd\ ^ \Bb\ and A Ki Ad¬ 
it is easy to show in this case that except near the 
edge of the disk (where edge effects change the values 
of and Sc becomes a constant multiple of the 
local disk eccentricity Cd- Indeed, if the precession of the 
disk is ignored, then Equation ( [2^ gives us that Cc = 
\Bd/Ad + ed\- Using Equations and § , this can be 
rewritten as 


— Od 


V’2 

2V-1 


1 


0.65ed, 


(30) 


as forp = 3/2, g = 1, we have V '2 ~ 1-82 and ipi « —0.55. 
In this regime, planetesimal survival is easier at larger 
semi-major axes simply because larger Op means lower 
ed{ap) for g > 0. 

4.3. Radial Behavior of 

The behavior of the relative planetesimal-gas eccen¬ 
tricity £iven by equations (251 and (29) is illustrated in 
Figure® which shows er as a function of orbital distance 
Op for different disk models. Unless otherwise noted, all 
panels assume the Kepler 16 binary parameters with our 
fiducial disk parameters as displayed in Panel A, and ap- 
sidal alignment of the binary and the disk, i.e. vJd = vjb- 
The most pronounced feature seen in all pane ls is the 
secular resonance, where A = 0, see Equation (26). It 


inevitably appears in the non-precessing disks at a a be¬ 
cause A changes sign there, as a result of At and Ad 
having different signs, see Figure The resonance lo¬ 
cation is different in panel B because of the lowered Eg 
in this particular disk model, which pushes oa out to 


4.8 AU, compared to a a ~ 1-9 AU in all other panels. 
Note that although formally Cc —>■ oo at the resonance, 
Cr stays hnite there. Moreover, the increase of at the 
resonance is often not very pronounced: unless dp > dc, 
the resonance has little effect on relative eccentricities 
because the planetesimals are closely coupled to the gas 
(see Equation (29)). The existence of the secular reso¬ 
nance previously suggested in the dra g-free environmen t 
by R13 was confirmed numerically by Meschiari (2014). 
In his simulations disk gravity excites the eccentricities 
of dp = 5 km planetesimals in a radially narrow disk re¬ 
gion around Op « 2.5 AU by a factor of ^ 10 compared 
to Cp at that radius in the absence of disk gravity. 

We see that lowering the disk eccentricity (panel C) 
causes to be dramatically lower outside of the secular 
resonance. This is because far from the star, disk gravity 
dominates, and we reach the limit of Section |4.2[ where 
Be is proportional to ed- On the other hand, reduction 
of eg has little effect near the star, where excitation is 
dominated by the binary. This is because very clo se t o 
the star, terms Ad and Bd are irrelevant in Equation (26). 
At the same time, disk eccentricity still enters through 
Cd, which is not negligible (see Equation (26)), so some 


variation of Op is still present near the star as eg is varied. 

Finally, panel D shows a shallower drop-off of planetes¬ 
imal eccentricity at large Op simply because we adopted a 
model with a different (lower) value of g, and Cp is in the 
disk-dominated regime outside the secular resonance. 


5. 


DYNAMICS WITH DISK-DOMINATED EXCITATION 
AND DISK PRECESSION 


We will find later in this work ((8.1) that it is gen¬ 


erally most promising to form planets outside of a few 
AU, relatively far from the binary. At large Op > a a 
planetesimal eccentricity excitation is dominated solely 
by the disk gravity, even though binary gravity may still 
affect planetesimal precession, see §3.1.1| 

For that reason we will now explore the limit in which 
the binary plays a negligible role in the planetesimal ec¬ 
centricity excitation, i.e. \Bb\ < \Bd\- At the same time, 
the binary is still allowed to contribute significantly to 
planetesimal precession, i.e. A = A/, + Ad in general. 
Thus, the current limit is analogous to the DD and DB 
regimes of SR15, with the addition of gas drag and disk 
precession. 

Here we also include a (realistic) possibility of the disk 
precession at some constant rate ibd- For the purposes 
of deriving an analytical solution, this rate must be inde¬ 
pendent of the semi-major axis, because the derivation 
of Rd in SR15 assumes fluid trajectories to be apsidally 
aligned at all Op. In practice this assumption is likely to 
break far from the binary. However, this is not a prob¬ 
lem as t he disk eccentricity is going to be very low there 
anyw ay (|Pelupessy & Portegies Zwart 2013 Meschiari 
2014|. 

Setting zud = Wdt and Hf, —)■ 0 in Equations ([l^-([l9l 


we find that, as shown in RS15a, these equations admit 
a periodic solution in the form 


rjBl 


-I 1/2 


COS (wdt -I- 4>) 

1 +{A-vJdyTj\ \sm{vJdt + (j)) 


(31) 


where the phase shift with res pect to the disk apsidal 
line is still given by equation (23) but with A replaced by 









































A — Wd- In this solution, the forced eccentricity vector 
rotates at the rate Wd and is fixed in the frame precessing 
with the disk. 


Solution (31) implies that the relative planetesimal-gas 
eccentricity is given by 


{A - Wd)Td 


— ‘-'C j - 1 

V 1 + [(^ - '^d)Tdf 

where the characteristic eccentricity Cc is now 

Bd 


— 


A - Wd 


+ Sd 


(32) 


(33) 


Also, equations ([27|)-(29) still hold with the provision 
that A is replaced witn A — Wd everywhere. 

5.1. Ejfects of Disk Precession 

We now explore the effect of disk precession on the 
behavior of planetesimal eccentricity. We examine the 
run of Cr in the disk for different disk models in Figure 
|5 assuming the binary parameters of Kepler 16. All 
disks have p = 1.5, q = 1, and Sq = 3,000 g cm“^. 
Although we show the behavior of starting at 1 AU, 
our disk-dominated excitation assumption is valid only 
for Op > Ob (the latter is shown by the green vertical 
line), which may be a problem for low cq (left panels in 
this figure) and has to be kept in mind. Because of that 
this figure would be accurate for all Op only if the binary 
eccentricity were zero. 

As a fiducial value of Wd we take the planetesimal pre¬ 
cession rate AjeAu at Op = 6 AU. This choice is almost 
entirely arbitrary, and is motivated only by the expecta¬ 
tion of the reduced disk eccentricity beyond this radius. 
Precession of the binary is likely to be an insignificant 
driver of the disk precession outside of a a, where the 
disk gravity dominates. For our fiducial system parame¬ 
ters, A|6AU = -1-2 X 10“^yr“k 

Figure 1^ shows that becomes independent of plan¬ 
etesimal size far from the binary, where Cc is low. This is 
due to apsidal alignment of planetesimals with the sizes 
shown in this Figure, as all of them are larger than the 
critical size in the outer disk, simply because disk ec¬ 
centricity is very low at large Op, see Equation (27). 


5.1.1. Secular resonances 

Equation ( [^ indicates that disk precession gives rise 
to two speciaTlocations in the disk. First, at the semi¬ 
major axis where 


A Wd^ 


(34) 


there is a secular resonance where Cc —> oo. This diver¬ 
gence happens because the relative precession between 
the planetesimal orbits and disk apsidal line vanishes, 
while the torque exerted by the non-axisymmetric com¬ 
ponent of the disk gravity is active. As a result, eccen¬ 
tricity can grow without bound in the absence of gas 
drag. This resonance is an obv ious generalization of the 
secular resonance discussed in §4.3[ 

Figure illustrates the non-trivial behavior of A, as 
it is a combination of Ad, and Aj,, which have opposite 
signs. We see that depending on the disk precession rate, 
there can be zero, one, or two secular resonances (34) 
associated with disk gravity. 



ttp, AU AU 


dp = .1km - dp = 1km - dp = 10km 


Fig. 4. — Run of relative planetesimal-gas eccentricity Cr in the 
disk dominated excitation regime, calculated for eight different disk 
models. We have assumed the fiducial system parameters except 
where noted. All models on the right have eo = 0.05, while all 
models on the left have lowered disk eccentricity, eo = 0.005. Dif¬ 
ferent panels correspond to different disk precession rates in units 
of AjgAU (also indicated in Figure^ — the planetesimal precession 
rate at 6 AU: = 0 (A,B), = -^Ioau (C,D), = ^leAU 

(E,F), ujd = 2A|6au (G,H). Depending on the disk precession di¬ 
rection and rate curves of e^ feature zero, one, or two secular res¬ 
onances. 

If Wd = 0 (non-precessing disk, yellow level in Figure 
Q , then there is only one resonance at the location where 
Wd -I- Ab = 0 (i.e. at 1.9 AU), and the situation is analo¬ 
gous to §4.3 see Figure |^,B. As in Figurethe larger 
planetesimals are excitea to higher eccentricities at the 
resonance because they are less damped by the gas drag 
(see Equation @).- 

Prograde precession gives rise to similar behavior, as 
shown in Figure [^,D. This is not surprising since Fig¬ 
ure [B shows that prograde precession (green level) sim¬ 
ply ^ifts the resonance location inwards, closer to the 
binary. This is indeed reflected in Figure |Jl, D, where 
the resonance is now closer to the star. 

As demonstrated in Figure retrograde precession 
can either remove the resonance if it is very rapid, 
\wd\ > |min(A)| (blue level), or give rise to two secular 
resonances as illustrated by the red level. Figure |^,H 
illustrates the former possibility. It is obvious that 
does not exhibit sharp features in this case. The mild 
bump around 3 AU is due to the reduced \A — Wd\ near 
the minimum of ACup). 

Finally, Figure |^,F shows the case of slower ret¬ 
rograde precession, with two conspicuous secular reso¬ 
nances at Qp ~ 2.2 AU and Op = 6 AU, in agreement 
with our expectations. It is clear that in this case 
can stay at a high level within an extended disk region 
between the two resonances, harming the prospects for 
planetesimal growth there. 

5.1.2. Regions of low Cr 
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A second type of special location in the disk is possi¬ 
ble when the disk processes in a prograde sense, i.e. in 
the direction opposite to the precession of the planetes- 
imal orbits. Equation (33) predicts that ec —>■ 0 when 


Binary precession rate 


Wd = A + Bd/cd- WhenThis condition is fulfilled, rela¬ 
tive velocities of particles of any size with respect to gas 
(and, consequently, also w.r.t. each other) vanish at this 
location, naturally promoting growth. This cannot oc¬ 
cur without disk precession, since Bd/cd is always larger 
in magnitude than A, except very close to the central 
binary where A has the same sign as Bd/cd- 
Using equations (|^ and §> , we see that Cc = 0 when 
Wd = Ah + Ad{l + '^/2'ipi), and outside a a the Ah term 
can be neglected. Since, by our assumption, thd is con¬ 
stant with radius, but Ad < 0 is not, one finds that in 
general a prograde disk precession {wd > 0) will give rise 
to a radius in the disk where Cc = 0, since ' 02 /V'i ~ 

This situation is clearly seen in Figure |^,D, where 
disk precession in the direction opposite to A creates a 
region around Op « 5 AU where —?> 0, independent of 
the disk eccentricity. This “valley of tranquility”, where 
relative planetesimal velocities are low, may represent 
a location where planetesimal growth is naturally pro¬ 
moted. 

To summarize this analysis, disk precession can be ei¬ 
ther helpful or harmful to planet formation depending on 
its direction and magnitude, and on the location in the 
disk. 

6. PLANETESIMAL DYNAMICS AROUND A PRECESSING 

BINARY 

One important difference between the dynamical en¬ 
vironments of P-type and S-type binary systems is that 
in the P-type systems, the precession rate of the binary 
itself can easily be comparable to or faster than t he plan¬ 
etesim al precession rates for much of the disk (Rafikov 
2013). It is therefore important to consider the preces¬ 


sion of the central binary in calculating planetesimal ec¬ 
centricities in situations where excitation due to the bi¬ 
nary is important. To that effect, we discuss diffe rent 
mechanisms driving binary precession in Section |6.1| and 
Appendix [B| Based on that, we then explore the role of 
binary precession using a simple disk model in Section 


6 .1. Binary Precession Rates 

We consider four drivers of binary precession: (1) 
general relativistic precession, (2) the quadrupole due 
to tidal interaction between the two stars, (3) the 
quadrupole due to stellar rotation, and (4) the gravity 
of the circumbinary disk. In Appendix ^ we provide 
estimates of the magnitude of each contribution. 

We find that precession due to the tidally induced 
quadrupole is dominant for massive stars that are close 
together, and disk-driven precession is dominant for less 
massive and more widely separated binaries. General rel¬ 
ativistic precession and rotationally induced quadrupole 
precession are always subdominant to one or the other 
of tidal or disk precession for the Kepler systems. 

Figure shows our estimate of the binary precession 
rate as a function of Mh and Ob- This assumes that 
fi = 1/3, t he a psidal mot ion constant k 2 = 0.13 (see 
Equations (B.2) and (B.3)) and the log o f the surface 
gravity (shown m the work of Claret (2012) to be nearly 


o 

Kepler 34 

n 

KIC 4862625 

• 

Kepler 16 

■ 

Kepler 413 

• 

Kepler 47 

■ 

Kepler 35 

u 

Kepler 38 

u 

KIC9632895 



I'l.yT 

Isxio^ 


i 


Fig. 5.— Binary precession rate (yr“^), as a function of aj and 
Mb- Thin white lines show the contours where disk precession 
equals tidal precession, and where it exceeds the latter by a factor 
of 10. Points correspond to known circumbinary systems listed in 
Table [l] We have assumed /i = 1/3 in calculating the precession 
rate. See text for more details. 


constant over a wide mass range for 1 Myr old stars) 
is 3.66 in cgs units. We have taken th ese numbers f rom 

'Claret (2012 | ) for 


the pre-main-sequence stellar models of < 
stars 1 Myr in age. We assume our fiducial disk model 
(Table in calculating the contribution of disk gravity 
to binary precession. 

We have additionally placed on this figure eight Kepler 
binary systems known to host circumbinary planets. The 
precession rates for these systems are only approximate, 
as they do not all have /i = 1/3 assumed in our calcu¬ 
lation, and there is a substantial change in stellar radii 
over the course of pre-main-sequence evolution. The con¬ 
tours correspond to locations in Ob — Mb space where disk 
precession rate equals tidal precession rate, and where it 
exceeds the latter by a factor of 10. If the disk cavity 
were larger (i.e. Oin were larger than 2ab), then the disk 
would drive substantially slower binary precession, see 
Appendix [B] 

We see from Figure that assuming precession to be 
disk-dominated is a good approximation for most of these 
systems, and particularly for Kepler 16, for which binary 
precession is dominated by the disk by a factor of hun¬ 
dreds over the tidal quadrupole precession. Throughout 
the rest of the paper, we assume that binary precession 
is solely due to the disk. 

Using Equations (1^ and (IB. 41), with -(/i = —0.55, we 
find 


^ =_ 0 . 15 ^ 

Ad Oin 


(35) 


for our disk model with p = 1.5. Equation (B.4) yields 
Wb because we are assuming ab/ai„ = 1/2, and 

Ad ~ a~^ from Equation ([^. Because Ad and Ab have 
different signs, |A| < \Ad\. As a result, even inside of the 
radius where Ad = thb, the major contribution to the 
relative precession between the binary and the planetes- 
imals may be coming from the binary precession — see 
Figure 
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6 .2. Axisymmetric Disk 

As an application, here we consider planetesimal dy¬ 
namics in the limiting case of an axisymmetric disk 
(eo = 0) with a central binary precessing at the rate 
Wh- T his section considers the same setup as in |Rafiko^ 
(|2013|), but with the addition of gas drag. 

Set ting Wb = Wbt and t 0, > 0 in Equations 


(18)-(19) we can easily solve them analytically. We find 
tnat tne solution for is identical to Equation (31) if we 
set Cd —?► 0 and replace Bd ^ Bb, xud —>■ xub in the latter. 
In this solution, the forced eccentricity vector rotates 
at the rate Wb and is fixed in the frame precessing with 
the binary. 

The relative planetesimal-gas eccentricity is still given 
by Equation (32) but with thd replaced by vjb- And 
instead of Equation (33) we now have 


I |. 

A — Wb 


(36) 


Given that binary precession is prograde, xub > 0, Fig¬ 
ure implies that there could be only one secular res¬ 
onance associated with binary precession (illustrated by 
the orange level calculated for our fiducial parameters of 
Kepler-16) and determined by the condition 


solution for the linear gas drag law described in Appendix 
[A] Examination of this solution shows that in general one 
^ould expect sec ular r eso nances of both types — given 


by Equations (34) and ( 37) — to exist in the disk. Given 
the discussion m §5.1.1 and |6.2| one expects up to three 
secular resonances to emerge, in particular, three reso¬ 
nances appear for slow retrograde disk precession, with 
two resonances being due to disk gravity (corresponding 
to the situation in FigureEE.F) and one due to the binary 
precession (see Equation p7[ )). This makes planetesimal 
dynamics even more complex than before. 

Neverthel ess, the gene ral features of the behavior 
outlined in §4.3| and |5.1| and shown Figures and re¬ 
main in place: eccentricity reaches high values at reso¬ 
nances, with the larger increase of e^ for bigger objects, 
less coupled to gas. At large separations disk gravity 
would still do mina te and dr ive to size-independent 
behavior, see §4.2| and |5.1.1| Lower disk eccentricity Cq 
would still result in lower planetesimal eccentricity, and 
so on. 

This completes our discussion of secular planetesimal 
dynamics in circumbinary disks. 

8. COLLISIONAL OUTCOMES AND GROWTH OF 
PLANETESIMALS 


A = 'djb- 


(37) 


It is located somewhat closer to the star than a a- Given 
all that, it is clear that the plot of Cr is going to be similar 
to Figure |^,D. However the “valley of tranquility” is 
going to be absent, and the secular resonance would now 
correspond to the location where the condition (37) is 
satisfied. 

In the case of very tight binaries (at, < O.lAU) binary 
precession is much faster than the planetesimal preces¬ 
sion rate A throughout the disk, thus strongly suppress¬ 
ing the exc itation due to the binary. This effect was first 
noticed in ( Rafikov|2013 ). Even more modest precession 
rates, such as those calculated for Kepler-16 in this pa¬ 
per (tbt, = 2.6 X 10“^yr“^) reduce the excitation effect 
of the binary by factors of a few. Looking at Figure 
we see that in most of the disk outside a a, the main con¬ 
tribution to the relative planetesimal-binary precession 
A — -ujb is due to the high value of Wb- Thus, binary 
precession effectively suppresses planetesimal eccentric¬ 
ity excitation due to its own non-Newtonian potential in 
the outer parts of the disk. 

7. DYNAMICS WITH BOTH BINARY AND DISK GRAVITY 
AND BOTH BINARY AND DISK PRECESSION 

To conclude our discussion of circumbinary planetes¬ 
imal dynamics we provide a general description of the 
Bp behavior in the most general situation when neither 
disk or binary gravity, nor disk or binary precession can 
be ignored, and quadratic gas drag damps planetesimal 
eccentricity. In this case the forced planetesimal eccen¬ 
tricity is no longer constant in time (even in some pre- 
cessing frame) and its amplitude varies. As shown in 
Beauge et al. (2010) and SR15, evolution of Bp can be 
viewed as a superposition of the two precessions in the 
eccentricity space, resulting in a limit cycle behavior with 
dp-dependent characteristics, which cannot be described 
analytically. 

Despite this complication, important insights into the 
general problem can be gained by analyzing the general 


To assess the prospects for planet formation in cir¬ 
cumbinary systems, we must couple our understanding 
of the planetesimal dynamics outlined in previous sec¬ 
tions with a prescription for the outcome of planetesimal 
collisions. We adopt that f rom RSI5b, who based thei r 
calculation on the results of Stewart & Leinhardt (2009). 
In their framework the outcome depends on the masses 
of the planetesimals, an assumption about their internal 
strength, and their collision velocity, which is determined 
by the Bp of each body involved in a collision. 

For completeness, details of our collisional prescription 
are reproduced in Appendix We categorize collisions 
in three groups based on the mass of the largest surviving 
fragment. 


• Catastrophic collisions leave a largest remnant con¬ 
taining no more than half the combined mass of the 
two colliding bodies. 

• An erosive collision leaves a largest remnant 
smaller than the larger of the two incoming plan¬ 
etesimals. 

• If the largest remnant is bigger than either of the 
two incoming bodies, we say that the collision leads 
to growth. 

Figure [^illustrates collisional outcomes in the space of 
sizes of colliding planetesimals di and ^2 at two different 
locations in the circumbinary disk of the Kepler 16 sys¬ 
tem, using the disk-dominated excitation approximation 
without precession, and with fiducial disk parameters. 
White contours enclose regions leading to catastrophic 
disruption. Black contours enclose regions leading to ero¬ 
sion. 

We see that collisions of bodies of exactly the same 
size lead to growth because collision velocities are small. 
Catastrophic disruption occurs only closer to the binary 
(panel A), since Cc is not high enough further out (panel 
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m/s 



d^, m 


Fig. 6. — Regions of catastrophic disruption and erosion for two 
different environments in the space of sizes of the two colliding 
partners. Panel A is made for our fiducial system at 3.5 AU with 
no disk or binary precession, assuming dynamics derived in ^ 
White lines enclose regions of destructive collisions for “strong 
(solid ) and “weak” (dashed) planetesimals (|Stewart & Leinhardt| 
|2009||. Black lines enclose regions of erosive collisions tor strong 
(solid) and weak (dashed) planetesimals. Panel B shows the same 
system at 6 AU. Here we see that although there is no destruction 
region for strong planetesimals, a large range of collisions still result 
in erosion. We have also shown the value of Sc and placed a white 
star at the locations of dc in each panel. Definition of the critical 
planetesimal size necessary for growth dmin is illustrated in the top 
panel. 

B) for catastrophic collisions to occur. We see that sim¬ 
ilar in size, but not exactly equal planetesimals tend to 
undergo catastrophic collisions, because their Bp are dif¬ 
ferent enough to result in significantly energetic colli¬ 
sions. However, even very unequal mass ratio collisions 
can lead to erosion. 

The white star in each plot corresponds to the size dc- 
We notice that this is near the region of catastrophic 
destruction in the plot. The lobes are biased towards 
dp > dc because dc is smaller than the size of plan¬ 
etesimal with the lowest critical velocity for destruction, 
which is about 100 m, see Figure 2 of RS15b. Because 
the destruction region is near dc, having a low value of dc 
means that planets can form starting from a population 
of smaller planetesimals. In fact. Figure [^suggests that 
if the planetesimal population were to contain only the 
objects with sizes > (30 — 50)dc, this population would 
have been completely immune to both catastrophic dis¬ 
ruption and erosion. 

In this paper, we consider catastrophic disruption to 
be the only obstacle to planetesimal growth. For a num¬ 


ber of reasons discussed further in Section |9.1[ we do 
not expect erosion to play a determining role in whether 
planetesimal growth occurs. 

In the following sections we examine the prospects for 
collisional growth in each of the dynamical regimes dis¬ 
cussed in Sections 1^-^ which provide different prescrip¬ 
tions for Cc and dc- Using these, we calculate the small¬ 
est planetesimal size dmin such that objects larger than 
dmin do not Suffer catastrophic disruption. This size is 
indicated in panel A of Figure This is the smallest 
planetesimal size that we can start from and grow larger 
objects without ever encountering catastrophic disrup¬ 
tion. Because of our neglect of erosion, overall planetes¬ 
imal growth requires only the existence of objects larger 
than dmin in the planetesimal populations. 

If this size is under 10 m, we conclude that planetesimal 
growth via collisional agglomeration is easy under those 
environmental conditions. The choice of 10 m is some¬ 
what arbitrary, but we note that changing from 10 to 100 
m in our plots would make very little difference to our 
conclusions. If dmin > 10 m, dmin provides an estimate of 
the minimum size of primordial planetesimals necessary 
to form planets in that environment. Unless otherwise 
specified, we are assuming the material properties ap¬ 
propriate fOT_fte “strong planetesimals” of [Stewart 


Leinhardt (2009). 


8.1. Planetesimal growth in the Disk-Dominated 
Excitation Regime 

We first consider the collision outcomes in the outer 
part of the disk where excitation from the binary is unim¬ 
portant compared with excitation from the disk, and use 
the results on planetesimal dynamics from Section 

Figure shows the size dmin, above which the growth 
is unimpraed by catastrophic disruption, as a function 
of semi-major axis and disk eccentricity for different as¬ 
sumptions about the disk precession. We have used our 
fiducial system parameters listed in Table Collision 
outcomes have been calculated using Equations ( [^ and 
(27) for Cc and dc- 

(Jne can see that in general, the outer region of the 
disk is favorable for growth, even starting from small 
(dp ^ 10 m) planetesimals. This is true for two reasons. 
First, local disk eccentricity is small at large Op, which 
leads to low values of Cc- Second, far from the binary 
dc oc Op ^ and becomes quite low (around a meter) at 
around 5 AU. This means that all the planetesimals have 
{A — Wd)Td ^ 1 and very near Cc- 


8.1.1. Effects of Disk Precession 

Looking at Figure [Tj we see that in the absence of 
disk precession, the outer region of the disk is friendly to 
planetesimal growth beyond 3-4 AU (for all planetesimal 
sizes down to 10 m) depending on cq. Starting with 
planetesimal sizes of a few km brings the inner edge of 
the growth-friendly region within 3 AU, even for disks 
with Co = 0.05 (twice the free-particle eccentricity). 

In the middle panel the growth region expands be¬ 
cause prograde disk precession (—AjeAu is positive) dra¬ 
matically lowers Cc- In this case, the disk precesses 
at vud = A\e,AU (A|6AU is the planetesimal precession 
rate at 6 AU), and we find nearly the whole outer disk 
(cLp > 3AU) to be conducive to planetesimal coagulation 
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Fig. 7.— Minimum planetesimal size dmin safe from catastrophic 
disruption in the disk dominated excitation regime (described in 
Calculation is done for our fiducial disk around Kepler 16, 
see Table Different panels correspond to different assumptions 
about the aisk precession rate and direction. The whited out region 
in the bottom left is where \Bf)\ > \B(i\ and we are not justified in 
using the disk dominated approximation. The dashed red line is 
drawn at the value of eo corresponding to the forced ecce ntricity 
of a free particle in the binary potential, see Equation ( |10[ |. 

even for highly eccentric disks. This is because of the 
valley of l ow e g in a prograde precessing disk, discussed 
in Section [5~T| and shown in Figure J^,D. 

On the other hand, retrograde diskprecession may give 
rise to a second secular resonance, which is very dam¬ 
aging to planetesimal growth, see Figure W- This res¬ 
onance makes conditions in the outer part of the disk, 
around 6 AU, much more hostile to planetesimal growth 
than in the absence of disk precession. In addition, com¬ 
paring panels A and C, we see that the inner resonance 
in panel C is also moved slightly outwards because of 
the disk precession. As a result, planetesimals within a 
broad radial interval of the disk (1-8 AU) end up being 
strongly dynamically excited. 


8.1.2. Effect of Planetesimal Strength 

Next, we consider changing the material properties of 
the p lane t esimals, i.e. changing the Qp n term in Equa¬ 
tion (|C.1|). iStewart & Leinhardt|(|2009| provide two pre- 

■ * j 10. I II • * I . *-1 


scriptions for depending on whether planetesimals 
are assumed to be solid rocks (strong planetesimals), or 


Fig. 8.— Same as Figure except t hat we are using the colli- 
sional prescription for weak aggregates JStewart fc Leinhardt|2009|l 
instead of that for solid rocks, see Appendix |tj| 

rubble piles (weak planetesimals). Figure shows the 
same dynamical environment as Figure but considers 
weak planetesimals. This does not make a big difference 
except for the sub-kilometer sized planetesimals, as there 
is not a big difference between weak and strong planetes¬ 
imals in the gravity-dominated regime (for dp ^ 1 km). 
Comparing Figures,!^ andwe see differences of several 
AU in the extent ofthe region where 10 m sized plan¬ 
etesimals are vulnerable to catastrophic disruption, but 
insignificant differences in the extent of the region where 
km-sized planetesimals are subject to destruction. 


8.1.3. Disk Mass 

Another model parameter which we vary is the surface 
density at 1 AU. Unlike the case of tight S-type systems, 
there is little reason to believe that circumbinary disks 
contain less mass than their counterparts around single 
stars. In fact, if anything, there seems to be e vidence 
for more m assive disks in circumbinary systems (Harris 

eTal lM^ . 

In Figure we consider both a denser and less dense 
disk than our fiducial system, with = 0, again using 
the dynamics discussed in Section We have adjusted 
the scale on the y-axis of Figure [Tin both subplots so 
that the bottom of the plot is near as, where excitation 
switches to being dominated by the disk, see Equation 
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Fig. 9. — Same as the top panel of Figure!^ but now for two 
different disk masses, resulting in different surface densities at 1 
AU indicated on panels. We use the system parameters in 
Table except for the disk density. 

( [I^ . Because qb depends on disk mass, the scales on 
the y-axes are not the same in the two panels. We have 
whited out the area in the plot where \Bb\ > \Bd\, as we 
do not have an analytic solution for Cc in that case. 

The main effect of changing the disk mass is the varia¬ 
tion in location of the secular resonance where A^ + Ai, = 
0: it moves out for lower Sg (and Md). For this reason 
lowering the disk mass is quite unfavorable for planet for¬ 
mation in the outer part of the disk. There may however 
be a region favorable to growth interior to the resonance, 
a possibility we explore further in Section 

For the very massive disk with Eg = 3 x 10“^ g cm“ 
we find that the dynamics look similar to the top panel 
of Figure (again accounting for the change in scale on 
the y-axisjT This is because in the disk dominated exci¬ 
tation regime, dc and Cc are i nde pendent of disk mass, 
see discussion after Equation (27). The only difference 
we expect from increasing the disk mass is that the secu¬ 
lar resonance moves inward. This is why the differences 
between Figures and become more striking as one 
moves inwards in semi-major axis. 

8 .2. CoUisional Outcomes in an Axisymmetric Disk 
with Binary Precession 


,-2 


As found in Rafikov (2013), a massive axisymmetric 
disk is very helpful tor reducing planetesimal collision ve¬ 
locities, both because of the enhanced planetesimal pre¬ 
cession and because of the induced binary precession. 
Not surprisingly, when including gas drag, this result re¬ 
mains valid, as we show now. 

In Figure we show dmin in tde axisymmetric disk 
approximation as a function of Eg and We are now 
using the dynamics discussed in Section]^ including the 
effects of non-z ero binary precession, in particular Equa¬ 
tions (35) and (36). We consider both strong (top) and 


Fig. 10.— Size dniim eis a function of Eg and ap for an axisym¬ 
metric disk with p = 1.5 in the Kepler 16 system. Planetesimal 
collisional velocities are calculated using the results of ij^ The two 
p anels correspond to the stron g and weak planetesimals discussed 
in|Stewart & Leinhardt|(|2009||. Binary preces sion rate ruf, is given 
as a function ot disk mass (Eg) by Equation (|B.4^. 


weak (bottom) planetesimals. 

Both panels exhibit similar structure, since the plan¬ 
etesimal dynamics are the same, independent of their 
material prop erti es. We see the secular resonance given 
by Equation (37) running diagonally across both pan¬ 
els. Its appearance is different from Figures [Tj and 
because those figures held disk mass constant. 

Exterior to the resonance, where A is disk- dom inated, 
dc is independent of disk mass, see Equation (27). How¬ 
ever a more massive axisymmetric disk lowers Cc by in¬ 
creasing the rate of relative precession between the bi¬ 
nary and planetesimal orbits, without adding to the ex¬ 
citation. We see that a high Eg allows unimpeded coag¬ 
ulation as close as 2.5 AU even for weak planetesimals. 
A low disk density is detrimental to planetesimal growth 
at large radii because the secular resonance moves out. 

Interior to the secular resonance, we see that a massive 
disk is actually harmful to the survival of smaller plan¬ 
etesimals. Indeed, fixing e.g. Op = I AU and increasing 
Eg leads to higher dmin- This is because in the inner disk 
A is dominated by the binary, A « A^, so that a more 
massive disk increases both the value of dc (because E^ 
increases in Equation (Hii), and Cc (by moving the reso¬ 
nance inwards). 

In the opposite limit of a low mass disk (Eg = 300 
g cm“^ corresponding to a disk mass of a few Jupiter 
masses), in situ planetesimal growth at Op « 1 AU is pos¬ 
sible with km-sized initial planetesimals. A light disk 
ensures that the secular resonance is at several AU, and 
therefore not playing a role in the dynamics inside of an 
AU. It also means that dc is still substantially smaller 
than a kilometer, even at an AU separation. For exam¬ 
ple, for an axisymmetric disk in our fiducial system with 
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Fig. 11.— Size (imin the inner part of the disk, adopting ap¬ 
proximation of no disk and binary precession which is v alid 

for small bodies. We calculate Cc by maximizing Equation l |26| l 
over mri — tuf,. White areas correspond to locations where Equa¬ 
tion ( |26[ | is no longer valid because > 1, for the largest size 

of planetesimal which is destroyed. System parameters are taken 
from Table [^except for the disk density. 

Eg = 300 g cm“^, the critical size is just 8 m at 1 AU. 

8.3. Planetesimal Growth Near the Star in the Limit of 
Short Stopping Times 

Planetesimal dynamics are complicated near the star 
because both gas drag and binary gravity are important, 
and the disk and binary do not process at the same rate, 
causing planetesimal eccentricities to be time dependent, 
see © Nevertheless, we still find analytic solutions for 
planetesimal orbits in two limits. One is the axisymmet- 
ric disk discussed in the previous section. The other is 
the limit of no disk and binary precession discussed in 
Section ID 

Planetesimals small enough to have stopping times 
short compared with the binary and disk precession times 
(max('d7b, tbdjrd ^ 1) should have their eccentricities ap¬ 
proximately described by Equations (20l-(24) with 
Wd given by the instantaneous orientation of the binary 
and the disk. We imagine \uJd\ to be smaller than \wh\, so 
we are only requiring W},Td ^ 1, when assessing the va¬ 
lidity of the approximation t hat stopping times are rapid. 

Planetesimal eccentricity (201 obtained in Section is 
a function of the mutual disk- binary apsidal orientation 
zud — zub- To be conservative in our estimate of the plan¬ 
etesimal destruction region, here we calculate the max¬ 
imum value of Cr as a function of zud — Wb in Equation 


(26), thus considering planetesimal dynamics in the least 


favorable part of the binary orbit. 

We use these assumptions to generate Figure which 
is similar to Figure We display as white the region 
where WbTd < 1 for planetesimals of size dmin- This is 
where we expect the approximation of a slowly-precessing 
binary to break down. 


We see a substantial difference of the outcomes depend¬ 
ing on the density of the disk. Because the gravitational 
perturbations are dominated by the binary (i.e. A is in¬ 
dependent of E(i), th e cr itical size dc is smaller in lighter 
disks, see Equation (l27|), leading to even km-sized bod¬ 
ies having more alignea orbits. Additionally, the secular 
resonance moves outwards for lower Eg, leading to lower 
Cc within an AU. 

As a result, in the case of a low-density, low- 
eccentricity disk, we see that it is possible to have plan¬ 
etesimals greater than a few km in size grow undisturbed 
by catastrophic disruption. This is a simpler scenario for 
planet esimal growth than the one described in |Meschiari| 
(20141 which relies on a pressure maximum to trap small 
dust to enable planetesimal growth. One caveat of the in- 
situ growth scenario with a low-Eg disk is that it may be 
difficult to grow Saturn-size circumbinary planets (such 
as Kepler 16b) because of the short supply of mass in 
such a disk. 

For the denser disks (higher Eg), catastrophic disrup¬ 
tion is inevitable unless nature creates i nitial planetesi¬ 


mals on the order of tens of km in size (Johansen et al. 


2012), see qO.2 These general conclusions are similar 
to those obtained in Section [8l2l 


Although our analytic results are formally accurate 
only for dp ^ several km, there is hardly a reason to 
believe that as they grow larger, planetesimals will be¬ 
come more vulnerable to destruction. Indeed, larger bod¬ 
ies are more resistant to higher velocity collisions. Also, 
we do not expect that they will have collision velocities 
dramatically higher than their lower-dp predecessors. In 
fact, it seems much more likely that large bodies will have 
smaller collision velocities. This is because close to the 
binary, |A| « \Ab\ 3> |d7f,|, so dc (for which \A\Td ~ 1, if 
we ignore in Equation (28l) is much smaller than the 
size of bodies for which TdWb ~ 1. Thus, larger bodies 
should be more resistant to destruction, see Figure]^ 


9. DISCUSSION 

Here we will revisit some of the assumptions that went 
into the model described in the previous sections, and 
explore what happens when they are relaxed. 


9.1. Erosion 

In the preceding sections, we have assumed that plan¬ 
etesimal growth is inevitable in the absence of catas¬ 
trophic disruption events. This is not necessarily the 
case, as frequent collisions with smaller objects can still 
lead to mass loss, even though none of them are severe 
enough to destroy the planetesimal. This is the erosion 
regime defined in SectionIn RSlSb (see their §4), we 
determined that in a collision between planetesimals of 
mass mi and m 2 , the critical velocity for erosion is in¬ 
dependent of m 2 in the limit that m 2 mi. Erosion 
can therefore be deleterious in regions where planetesi¬ 
mals are safe from catastrophic disruption because a low 
value of dc does not have as large of a protective influ¬ 
ence. 

If characteristic eccentricity Cc is high enough for col¬ 
lisions to be erosive in the limit of vanishing collision 
partner mass, then erosion might be expected to inhibit 
planetesimal growth. Looking at Figure we see that 
erosion is apparently ubiquitous, even in a regime where 
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strong planetesimals do not suffer catastrophic disrup¬ 
tion. Therefore, if a substantial fraction of the mass that 
a planetesimal encounters is in small bodies, it might 
have difficulty growing. This outcome could be avoided 
if the majority of the mass of solid material in the disk 
is contained in objects with sizes larger than dmin) or if 
collision rates between large and small bodies are sup¬ 
pressed. 

In practice, we do not expect erosion to be the major 
obstacle to planetesimal growth because small bodies are 
likely to be rapidly flushed out of the system due to gas 
drag in the slightly sub-Keple rian disk. 


Indeed, using the results of Adachi et al. (19761, and 


ignoring the relative particle-gas eccentricity, we estimate 
an in-spiral timescale of 


256 \/^ iPpdp 


dp 507^^(73/2)'"^ 

« 6 X lO^yrx 

Pp dp 3,000 g cm“^ 

3 g cm“3 1.3m Eq 




(38) 


Mb 


0.89Mq 


1.5 


9/4 
^p,5 ■ 


In this estimate we have taken for dp the critical size of 


1.3 m given by Equation (27) at 5 AU. Particles 


with sizes above dc are less likely to erode larger bodies 
due to their aligned orbits. This is a lower bound, as 
non-zero gas-particle eccentricity will only accelerate in¬ 
spiral. 


Furthermore, as discussed in Weidenschilling (1977), 
for these small bodies moving slowly relative to the gas, 
the drag is actually not in a quadratic regime, and is 
stronger than predicted by our assumed quadratic drag 
law, making the in-spiral more rapid. The particles 
which in-spiral most rapidly are those for which UpTd « 1. 
Using our quadratic drag law, ripTd = 1 at 5 AU for parti¬ 
cles with sizes of about a centimeter. The fact that we are 
doing this estimate at 5 AU, instead of the conventional 
1 AU, leads to meter sized particles lasting substantially 
longer than t he 100 yr timescale ge nerally discussed in 
the literature ( Weidenschilling||l977). 


The significance ol the estimate (38) is that the disk 
cannot store most of its mass in bodies smaller than dc 
because Tm is considerably shorter than the several Myr 
disk lifetime. As a result, the disk would rapidly lose all 
of its solid material by radial drift towards the central 
binary. By assuming most of the mass in the disk to be 
in planetesimals larger than several meters initially, we 
thus make erosion by smaller bodies to be a subdominant 
effect. 

In addition, very small bodies with stopping times 
shorter than the dynamical time n~^ orbit with sub- 
Keplerian velocities due to their strong coupling to gas, 
which experiences radial pressure support. These objects 
would impact larger bodies at speeds of several tens of m 
s“^, even in the absence of any secular excitation. The 
same would happen in protoplanetary disks around sin¬ 
gle stars, which according to recent statistics of exoplan¬ 
ets, have no diffuculty of forming planets. Thus, erosion 
by very small objects can be overcome in circumbinary 
systems in the same way as this happens around single 
stars. 


9.2. Growth via “Lucky” Bodies 

It is likely that planetesimal growth can occur even 
in presence of some catastrophic disruption, since the 
detrimental effect of these can be offset by other favor¬ 
able collisions which enable growth. This balance of the 
mass loss and gain has be en explored for th e grow th of 
centimeter-sized bod ies by Windmark et al. (2012) and 
Garaud et al. (2013). These authors consider a distribu¬ 
tion ot encounter velocities and collision outcomes, rather 
than assuming all encounters to be at the mean velocity. 
By using such statistical approach they find over an order 
of magnitude change in the size of the largest particles 
that can grow in their coagulation simulations. Thus, 
statistical nature of coagulation may be an important 
factor of planetesimal growth, which we did not account 
for here. 

In the km-sized planetesimal regime, bodies generally 
become both more resistant to collision, and experience 
lower velocity collisions with similar sized objects as they 
grow large r bec ause their size moves further from dc^ see 
Equation (C.4). It therefore seems likely that including 
the full distribution of collision velocities will make an 
even larger impact in this scenario than it did for the 
growth of centimeter sized grains, since a km-sized body 
that grows via a few lucky collisions becomes harder to 
destroy. We leave the detailed exploration of the statis¬ 
tical growth of planetesimals in circumbinary systems to 
future work. 

9.3. Limits of the Apsidally Aligned Regime 

In the previous sections we have assumed that free ec¬ 
centricity of planetesimals has been completely damped 
by gas drag. As the free eccentricity goes away on the 
timescale r^, we need to demonstrate that is shorter 
that the characteristic planet formation timescale. 

At a gi ven semi-major axis in the disk, we can use 
Equation (17) to solve for the critical damping size ddamp 
at which the free eccentricity damping time is equal 
to some characteristic time r (not to be confused with 
the critical size dc for which the damping time is roughly 
the precession time A“^!). If we set r equal to the ex¬ 
pected lifetime of the circumbinary disks (Myrs), then all 
objects with dp > ddamp will not have a chance to settle 
to their quasi-equilibrium forced eccentricities during the 
disk evolution. This woul d vio late our basic assumption 
of damped Cfree stated in §3.4| and would make planetes¬ 
imal growth more complicated. 

To demonstrate this last po int, we repeat the calcu¬ 
lation illustrated in Figure 10 however, now we do not 
assume planetesimal orbits to be aligned according to 
their forced eccentricity values. Ins tead, we set ei 2 to be 
given by Cc instead of by Equation (C.3). In the approx¬ 
imation where planetesimals decouple instantly from the 
gas with eccentricity equal to e^, their eccentricity vector 
Bp circulates around the forced eccentricity with a mag¬ 
nitude |ed - eforcedi- As shown in RS15a, |ed - eforcedl 
is just Cc. Therefore a typical relative eccentricity ei 2 
between planetesimals in this population will be of orde r 
Cc- This is the approximation explored in Rafikov ( 2013[). 

Results of such calculation are shown in h'igure in 
Comparing Figures [T0| and [T^ shows that the alignment 
assumption makes a substantial difference, with the zone 
of catastrophic disruption significantly extended in Fig- 
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Fig. 12.— c?min for tho same dynamical environment as Figure fTo] 
but assuming planetesimal orbits to be unaligned so that ei 2 = ec- 
Comparison of these two figures shows that the apsidal alignment 
is a very important effect, greatly facilitating planetesimal growth. 


ure |12| compared with Figure [T0| for a given planetesimal 
size. Thus, lack of free eccentricity damping resulting in 
apsidal misalignment of planetesimal s should have a dele¬ 


teriou s effect on planetesimal growth (Bromley & Kenyon 


20151. 


We expect protop lanetary disks to last for a few Myr 


(Haisch et al. 


2001). This motivates us to set a charac¬ 


teristic timescale in the deter min ation of the critical size 
ddamp to be ^ Myr. Figure 13 shows ddamp for r = 3 
Myr as a function of Op for several sets of disk parame¬ 
ters. This size is calculated numerically from Equations 
0 , (0, (0, pl), g, p]), and 
^We see that genei^y witninK) ATfothe damping time 
is shorter than the disk lifetime for kilometer sized plan¬ 
etesimals, depending, of course, on the density of the 
disk. Although outside of 5 AU, we may be concerned 
about the alignment of 10 km planetesimals, at these 
radii, such large planetesimals are generally well above 
the mass threshold for collisional growth (see |6.2[ ). This 
suggests that we are justified in considering the planetes- 
imal orbits to be aligned. It is also worth noting that a 
higher value of Cc leads to faster apsidal alignment, so 
the environments with high Cc, which present the highest 
danger to planetesimal coagulation are the same environ¬ 
ments which are helped most by the apsidal alignment. 

Collisions between planetesimals may also cause mis¬ 
alignment. However in the regions of the disk with long 
alignment timescales, planetesimals greater than ^ 10 m 
in size possess nearly apsidally aligned equilibrium or¬ 
bits, so collisions between them are unlikely to lead to 
much misalignment. 


9.4. Disk Density Structure 

One may wonder how our results would change if the 
slopes of the Ed(ap) and Cd^ap) dependencies were var¬ 


Fig. 13.— Maximum size of planetesimal in th e d isk dominated 
excitation regime (where Cc is given by Equation ( |30| l) such that 
is less than 3 Myr. We have assumed fiducial system parameters 
except where noted. 


ied. Observational evidence for the values of p and q is 
rather scant and is based predominantly in the sub-mm 
dust continuum observations on scales (tens of AU) much 
larger than the ones considered in this worl|^ (Andrews 


et al.|[2009). However, the total masses of the circumbi- 


nary di sks are found to be in the range of several percent 
of MoJlHarris et al. 2012), in agreement with our esti¬ 
mate (01 

We nave picked our fiducial value of Sq to match an 
overall disk mass of a few percent of Mq, believing this 
to be better constrained than the surface density Eq at 
1 AU. If we change p while keeping the overall mass of 
the disk constant, we see that Eq would vary with p as 


So(x(2-p)(^^y. (39) 

Equations 0 and 0 show that E^; and '0i are the only 
disk properties thar determine Ad- The coefficients ipi 
and i /’2 in Equations 0 and 0 change by less than a 
factor of 2 between p = 1.5 and p = 0.5, so we see that 
at a given radius, changing p at constant disk mass has 
about t he sa me effect as changing Eq, as was done in 
Section [8.1.3[ In other words, because ipi and i /)2 vary so 
weakly with p and g, what matters most for the disk dis¬ 
turbing function is the local value of the surface density, 
not the distribution in the whole disk. And it was shown 
in SR15 that the edges of the disk are unimportant in 
determining the disturbing function for the values of p 
and q that we consider in this paper. 


9.5. Non-Secular Terms in the Disturbing Function 


® ALMA may be close to marginally resolving these scales in the 
future. 
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In this study, we ignore the possibility of mean-motion 
resonances between the bin ary and the planetesimals. 
Previously, Meschiari (2014) found in simulations of the 
Kepler 16 system that planetesimals got trapped at the 
5:1 resonan c e. Sim ulating the S-type 7 Cephei system, 
Leiva et al. (2013) find planetesimals getting trapped in 
hrst order resonances as high as 16:1. While it is not 
clear if these results would hold far from the binary where 
gravity of a massive eccentric disk dominates, the mean 
motion resonances may still play an important role in the 
dynamics inside of 1 AU. 


Op < 1 AU Meschiari (2014) finds that km-sized plan- 


It was also suggested in |Paardekooper et al. (2012) 
that short-period terms in the disturbing function, vary¬ 
ing with the binary period or the planetesimal orbital pe¬ 
riod will add additional eccentricity. These seem unlikely 
to mis-align aligned orbits: short period terms should af¬ 
fect planetesimals of all sizes equally, not changing their 
relative velocities, as they act on a much shorter time- 
scale (~ n~^) than the gas drag damping time for 
the sizes we are concerned with. It was furthermore 
shown both an alytically and in simulation by [Bromley fc| 
Kenyon (2015) that in the Kepler 16 system these terms 
only induce absolute eccentricity variations of about a 
percent at 0.7 AU (relative variation should be much 
smaller). 

10. COMPARISON WITH PREVIOUS WORK ON PLANET 
FORMATION IN BINARIES 

Here we compare our findings with existing results on 
planet formation in both circumbinary and circumstellar 
systems and put them in broader context. 

10 . 1 . Previous Work on Circumbinary Planet 
Formation 

Despite the different setup and the range of physical ef¬ 
fects taken into account in this work, many of our conclu¬ 
sions are similar to those reached in a number of previous 
st udies of planet formation in circ umbinary systems . 


Paardekooper et al. (2012) and Meschiari (2012) nu- 
merically explored the interactions between swarms of 
planetesimals embedded in an axisymmetric gas disk 
around a central binary. They included the gravitational 
perturbations from the binary, and gas drag from the 
disk, but not the gravity of the disk, which as we know 
now (Rafikov 2013, SR15, RS15a) provides the most im¬ 
portant gravitational effect outside of a few AU. They 
both found it difficult to explain in-situ accretion at 
< 1 AU separation without invoking initial planetesimals 
g reater than 10 km in size. 


Meschiari (2012) concluded that formation outside of 
4 AU was possible starting from km-sized planetesimals, 
which is only slightly less optimistic than our conclusions 
for the axisymmetric disk of a similar mass. Our inclu¬ 
sion of disk gravity in this paper pushes the boundary 
of the coagulation region only a little bit inward because 
we find that the secular resonance, emerging when we 
account for disk gravity, makes conditions in the disk 
a round 2 AU unfav orable for planetesimal growth. 


Meschiari (2014) did include disk gravity in some of 
his simulations and found prominent secular resonance 


predicted in Rafikov (2013). However, he was primar¬ 
ily interested in in-situ growth of circumbinary planets 
within 1 AU and did not explore planetesima l gro wth 
at several AU in as much detail as we do in 98.1| At 


etesimals can grow by accretion of small collisional de- 
bris. By carefully examining planetesimal dynamics in 
§8.3| we find that growth starting with km-sized objects 
is possible even with more standard collisional scenario 
not involving accretion of small particles, as long as the 
disk density is not very high and the characteristic plan¬ 
etesimal size dr is small. 


Simulations of Marzari et al. (2013) also include the 
gravitational effect of a non-axisymmetric disk. They 
calculated the disk structure in 2D geometry, and found 
complex behavior for disk surface density and eccentric¬ 
ity with radius: their disk maintains substantial Cd (sev¬ 
eral percent) out to nearly 10 AU. This eccentricity pro¬ 
file looks very different (much higher) than found in simu¬ 
lations of Meschiari (2014), Pelupessy & Portegies Zwart 
( 2013), and may reile ct transient behavior. As a result, 
Marzari et al. (2013) found it difficult to grow even 25 
km planetesinials out to 10 AU because of eccentricity 
excitation due to disk gravity. This is at odds with our 
findings simply because we think that such eccentricity 
profiles are unlikely and assume to decay with radius. 
Another reason for the discrepancy is the duration of 
their runs (~ 10 "^ yr), which is likely shorter that Td for 
the 5-25 km planetesimals they are considering. This 
leaves significant undamped free eccentricity and results 
in large collisional velocities even among planetesimals of 
the same size. 


Bromley & Kenyon (2015) performed an analytic study 
01 planetesimal dynamics, and showed in particular that 
short-period terms in the disturbing function do not re¬ 
sult in increased collision velocit ies of planetesim als, in 
agreement with the argument of Rafikov (2013). How¬ 
ever, their prime focus was on pointing out the existence 
of a set of non-crossing, aligned orbits of planetesimals, 
which they found explicitly neglecting gas drag. They ar¬ 
gued that such orbits would allow planetesimals to collide 
with very low relative speeds (as low as in axisymmetric 
disks around single stars) and grow efficiently. 

Based on our results, we interpret these trajectories as 
orbits with fully damped free eccentricity and only forced 
eccentricity remaining, which by itself requires some ef¬ 
fective damping process (most likely gas drag) and w orks 


only for objects with dp < 10 km at several AU, see (9.3 
For such objects gas drag would in fact make forceo^; 
size-dependent (0 leading to substantial collisional 
velocities for planetesimals of different sizes even when 
all free eccentricity is damped, as clearly demonstrated 
in our (e.g. see Figure]^. This dynamical size seg¬ 
regation was previously broadly discussed in t he context 
of planetesim al growth in S-type systems (e.g. Thebault 
et al. (2008), RS15a). Additionally, when one consid¬ 
ers disk gravity, secular resonances arise, leading to orbit 
crossing even assuming apsidally aligned orbits, due to 
the rapid radial variation of Cp. 

10 . 2 . Large Initial Planetesimals 

It is possible that some mechanism other than colli¬ 
sional agglomeration produces a population of large ini¬ 
tial planetesimals resistant to collisional d estruction by 
virtue of the i r size. It has been proposed (Goodman & 
Pindor 2000 Youdin & Goodman 2005) that coupling 
between solid material and the disk could produce an 
instability that would lead to overdense rings of mate- 













































































18 


rial which would collapse to form planetesimals. This 
so-calle d streaming; instability has been explored numer¬ 
ically in Johansen et al. (20121, and was found to lead to 
rapid formation ot planetesimals several hundreds of km 
in size. 

The size distribution of minor bodies in our Solar Sys¬ 
tem provides clues as to the size of the initial planetesi¬ 
mals, however studies of this have come to conflicting re- 
(|2009|) find evidence for 100 -1000 
Is in the current size distribution of 


suits. Morbidelli et al. 
km initial planetesima. 


asteroids. However, |Weidenschilling (|2011|), using a dif¬ 
ferent accretion model, is able to reproduce the current 
size distribut ion starting from <0.1 k m planetesimals. 
Additionally, Schlichting et al. (20131 studied the size 
distribution of objects in the Kuiper belt, and came to 
the conclusion that initial planetesimals on the order of 
a km in size provide the best fit. 

Our present work shows that at least in circumbinary 
systems, rapid formation of large {dp ~ 10^ km) planetes¬ 
imals is not necessary and planetesimal growth is possible 
starting from km-sized (or even smaller, at large separa¬ 
tions) bodies. Thus, at present, circumbinary exoplanets 
do not provide a strong argument in favor of stream¬ 
ing instability being the dominant planetesimal forma¬ 
tion mechanism. It is also unclear whether the streaming 
instability would function in the perturbed circumbinary 
environment. 


10.3. Differences Between P-type and S-type Binaries 

This paper has a lot in common with RS15a,b. In this 
section we highlight some differences of the physics that 
goes into understanding dynamics in P-t ype systems. 

In S-type systems, there is evidence (Muller & Kley 
20121 that massive disks end up apsidaliy aligned with 
the binary star. This alignment gives rise to a dynami¬ 
cally cold region in the protoplanetary disks of the S-type 
systems discovered in RS15a, which cannot exist in the 
circumbinary configuration. Circumbinary disks can in 
fact contain a dynamically cold region, but only if they 
are precessing in the prograde sense (see panels C and D 
of Figure]^. 

The precession of the central binary cannot be ignored 
in P-type systems. This helps to lower the planetesimal 
eccentricity excitation due to the binary. This effect is 
absent in S-type systems. 

Circumbinary disks can be much more extended than 
circumprimary ones, simply because they are not tidally 
truncated by the companion on the outside. For that 
reason planet formation can (and is likely to) occur at 
large radii (with subsequent inward inigration of grown 
planets). We see from Equation ( p7| ) that everything 
else being equal, the critical size dc at which collision 
velocities between similar sized planetesimals are highest 
scales as dc Op The low value of dc in the planet¬ 
forming region leads to the conclusion that all bodies 
km-sized or larger are very nearly apsid aliy aligned as 
dp S> dc for all of them, see Equation (C.4). Additionally, 
the Keplerian velocity is smaller at higher separation, 
leading to further reduction of collision velocities. This 
permits planetesimal coagulation starting from relatively 
small initial sizes, 10 - 100 m, even in regions of higher Cc 
than is permitted in the S-type systems where massive 
planets reside at 1-2.5 AU. 


Another consequence of the low density environment 
at large separation in circumbinary systems is that for 
a substantial fraction of the likely planet-forming part 
of the disk, there is some concern that alignment with 
the gas may not occur within the disk lifetime, see Sec¬ 
tion |9.3[ This was not an issue for the S-type systems 
because their compact protoplanetary disks were likely 
much denser around 1-2 AU, effectively damping free ec¬ 
centricity and aligning planetesimal orbits. 

11. SUMMARY 

We studied planetesimal dynamics in circumbinary 
disks with the goal of understanding the conditions lead¬ 
ing to planetesimal growth, which is a natural step to¬ 
wards formation of circumbinary planets such as Kepler- 
16. Our study simultaneously considers (1) gas drag (in¬ 
cluding non-trivial radial pressure support in the gaseous 
disk), (2) gravity from the eccentric precessing binary, 
and (3) gravity from the eccentric precessing disk. We 
found analytical solutions for planetesimal eccentricity 
behavior in many important limits. 

We estimate the precession rate of the central binary 
and find binary precession to play a non-trivial role in 
the determination of planetesimal dynamics. We believe 
erosion by small bodies to not be a major hindrance to 
growth, as such objects should spiral in towards the cen¬ 
tral binary on timecales of tens of thousands of years. 
Based on our analytical results, we make the following 
conclusions: 

• We find disk gravity to play the dominant role in 
planetesimal dynamics outside of several AU. Sec¬ 
ular resonances (up to three for some choices of 
the disk and binary precession rates) significantly 
complicate planetesimal dynamics at Op ^ several 
AU. 

• Apsidal precession of the central binary reduces the 
direct effect of binary gravity on planetesimal or¬ 
bits. For many circumbinary systems discovered 
by Kepler precession is dominated by the gravity 
of the protoplanetary disk. 

• If planet formation is precluded only by catas¬ 
trophic disruption, then we find that forming 
Kepler-16 and similar planets in situ requires that 
initial planetesimal sizes were large and disk masses 
were small. For example, even in the most favor¬ 
able case of an axisymmetric disk with a mass of 
only a few Jupiter masses, we still require 1 km ini¬ 
tial planetesimals to avoid catastrophic disruption 
at 0.7 AU. If the disk were eccentric or more mas¬ 
sive, then even larger initial planetesimals would 
be required. 

• Formation outside of ^ 3 AU is much easier and 
more likely. Here a dense (E(i(lAU) > 10“^ g 
cm“^) disk (in agreement with sub-mm observa¬ 
tions) is helpful as its gravity moves the secular 
resonance inwards of 2 AU, lowering planetesimal- 
planetesimal collision velocity. The dense disk also 
provides enough eccentricity damping to align ~ 10 
km-sized planetesimals out to 10 AU during the 
disk lifetime. 
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• Disk precession can either facilitate or hinder 
plantesimal coagulation, depending on the direc¬ 
tion and magnitude of precession. Slow retrograde 
precession results in emergence of a destructive sec¬ 
ular resonance at several AU. On the other hand, 
prograde precession leads to a dynamically favor¬ 
able location in the disk where the forced eccentric¬ 
ity is the same as the local gas eccentricity, relative 


planetesimal eccentricities are small, and there is 
no dynamical barrier to coagulation. 

This work thus provides a dynamically-motivated pic¬ 
ture of planetesimal growth towards building planetary 
cores in circumbinary protoplanetary disks. 

This work is supported by NSF grant AST-1409524. 
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APPENDIX 


PLANETESIMAL DYNAMICS WITH LINEAR DRAG 


In the case of linear gas drag (i.e. independent of e^) we can obtain the complete analytical solution for planetesimal 
dynamics including the precession of both the binary and the disk. With fixed and Wbit) = vubt, Wd{t) = Wdt our 
general evolution equations (18)-(19) admit the following analytical solution: 


f-freef- 


-t/Td 


J COS {At + Wq) 
[sin {At -I- wq) 



(A.1) 


Here Cfree and zuq are the free eccentricity and periastron angle, ^ and hf^d are the components of the forced 
eccentricity associated with the disk, and kf^b-, and hf^b are the components associated with the binary. These are 
given by 


and 



gg+^d^d 1 / cos {7JJd{t) + (jid) \ 
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Fig. 14.— Dependence of the dimensionless factor p.) in the disk-driven binary precession rate d7disk (see Equation l |B.4[ |) on 

the relative size of the inner cavity ain/a^ and binary mass ratio /i. 


The relative particle-gas eccentricity is given by 


kf,b 

^f,b 


^ Bd+eg{A-dod) f COS {wd{t) — (j)^) 

[sin {zUd{t) - (pr) 


COS0r = 


TdiA-T^d) 


[l-rr|(^-d7d)2] 


1/2 • 


(A.4) 


It is clear that is in general a function of time. 

BINARY PRECESSION RATE 

Here we summarize results on the four major causes of binary precession which were discussed in Section We 
write the precession rate as tbf, = 'djoR + + '^rj) + i^diski where wqr, vjtj: '^rjj djdisk are the precession 

rates due to general relativity, tidal and rotational stellar quadrupoles, and disk gravity, respectively. 

General Relativistic Precession 

To first order in eccentricity, precession of Keplerian orbits due to general relativity is (Misner et al.|[l9^) 


'^GR = 


HGMt) 

c^al'^ 


1.5 


= 6.9 X 10"V"^ 




0.89Mo 
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0.22Ht7\ 
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(B.l) 


Precession Due to Stellar Quadrupoles Induced by Tides and Rotation 


P recession r ates d ue to quadrupoles induced by tidal forces and from the rotational bulge are given in Sterne (19391 
and Shakura (1985). They are given by 


'^T,j = I5k2jnb 
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Precession Due to Disk Gravity 


To calculate the binary precession due to disk gravity, we use Equations (20) and (A3) from Rafikov (2013), which 
assume that power law dependence of the disk surface density is sharply truncated at the inner radius Oin. These state 


"^disk — M)^6 


SoOp^b 

M^a. 


i+p 


o 1 / 0.89iVfQ 
2.6 X 10"^yr-i ——^ 

V ^b 


0.5 

f—) 

V0.22/ 


1/2-p 


3,000 g cm “2 0.46 


(B.4) 


Here 4>{flb/o,in) is a function of the binary mass ratio /i and ratio of binary semi-major axis to the radius of the central 
gap in the disk ab/a-m, and we have assumed ab/a-m — 0.5. Figure 14 shows that ^ is a slowly varying function of 
Qb/ain, so roughly zbdisk ~ {a-b/■ For p = 1.5, assuming ab/oin to be 1/3 instead of 1/2 would lower Wb by 
almost a factor of 3. 


COLLISIONAL OUTCOMES 


Consider two bodies of mass mi and m 2 , colliding at speed Ucoii- Stewart & Leinhardt (2009) give the mass of the 
largest remnant Mir as 


Mir 

Mtot 


-0.5(Qfl/Q*fli5 - 1) + 0.5 


(C.l) 


Here, Mtot = Wi -I-m 2 is the total mass of the colliding bodies, Qr = 0.5mim2u/oii/Mtot center of mass specific 

kinetic energy, and Q*rjj is a quantity dependent on the material properties and Mtot- As shown in Figure 2 of RS15b, 
the critical velocities for catastrophic disruption are on the order of several m s“^ for dp = 100 m planetesimals. The 
critical speed reflects both the strength of the body, and reaccumulation of fragments due to gravity. Material strength 
becomes subdominant to gravity for bodies larger than a few hundred meters. 

We take collision velocity to be given by 


^coll 




2G(mi -I- m 2 ) 
di (^2 


(C.2) 


where the first term is due to the relative velocity at infinity, 612 = \ep{di) — ep(d 2 )|, and the second is due to the 
potential energy of the colliding objects. It should be noted (RS15a) that the actu al velocity at infinity between two 
colliding bodies can be anywhere between 0.5ei2VK and ei 2 VK- Thus Equation (C.2) may overestimate the true Ucoii by 


up to a factor of 2. Because the collision velocity must be several times the escape velocity in order to be destructive, 
the potential energy term is relatively unimportant. 

We can express the relative eccentricity of two planetesimals in terms of Cc, dc and the planetesimal sizes di and ^2 
as (RS15a) 

„ |(A-d7d)(ri-Ta)! 


' y^{l + {A- Wd)‘^Tf){l + (A - ti7d)2r|) 


Here, |(A — Wd)Ti\ {i = 1, 2) is given in terms of dc and di by Equation (28). 

In the commonly encountered limit of di and ^2 both being much greater than dc 


this simplifies to 


ei2 


Ml - d2\dc 
did2 


(C.4) 


We see from this that for encounters between two large bodies, the relative eccentricity is never larger than 
ecdc/niin(di,d2). 




















